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Abstract 

In  this  paper  we  explicitly  construct  the  entropy  solutions  for  the  Lighthill-Whitham- 
Richards  (LWR)  traffic  flow  model  with  a  flow-density  relationship  which  is  piecewise  quadratic, 
concave,  but  not  continuous  at  the  junction  points  where  two  quadratic  polynomials  meet, 
and  with  piecewise  linear  initial  condition  and  piecewise  constant  boundary  conditions.  The 
existence  and  uniqueness  of  entropy  solutions  for  such  conservation  laws  with  discontinuous 
fluxes  are  not  known  mathematically.  We  have  used  the  approach  of  explicitly  construct¬ 
ing  the  entropy  solutions  to  a  sequence  of  approximate  problems  in  which  the  flow-density 
relationship  is  continuous  but  tends  to  the  discontinuous  flux  when  a  small  parameter  in 
this  sequence  tends  to  zero.  The  limit  of  the  entropy  solutions  for  this  sequence  is  explicitly 
constructed  and  is  considered  to  be  the  entropy  solution  associated  with  the  discontinuous 
flux.  We  apply  this  entropy  solution  construction  procedure  to  solve  three  representative 
traffic  flow  cases,  compare  them  with  numerical  solutions  obtained  by  a  high  order  weighted 
essentially  non-oscillatory  (WENO)  scheme,  and  discuss  the  results  from  traffic  flow  per¬ 
spectives. 
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1  Introduction 


Lighthill  and  Whitham  (1955)  and  Richards  (1956)  independently  proposed  a  macroscopic 
model  of  traffic  flow  to  describe  the  dynamic  characteristics  of  traffic  on  a  homogeneous  and 
unidirectional  highway,  which  is  now  known  as  the  LWR  model  in  the  literature  of  traffic 
flow  theory.  Although  a  substantial  amount  of  work  has  been  conducted  to  improve  the 
modeling  approach  of  traffic  flows  in  many  directions,  the  LWR  model  is  still  widely  used 
for  the  modeling  of  traffic  flow,  because  of  its  simplicity  and  good  explanatory  power  to 
understand  the  qualitative  behavior  of  road  traffic.  The  results  that  are  obtained  from  the 
LWR  model  are  generally  adequate  for  many  applications  such  as  traffic  management  and 
control  problems. 

The  LWR  model  is  formulated  as  a  scalar  hyperbolic  conservation  law  and  is  often  solved 
by  finite  difference  methods  (Daganzo,  1995;  LeVeque,  1992;  Lebacque,  1996;  Michalopoulos 
et  ah,  1984;  Wong  and  Wong,  2002a;  Zhang  et  al.,  2003).  The  main  difficulty  in  designing 
efficient  and  high  order  finite  difference  methods  for  the  LWR  model  or  in  general  for  hy¬ 
perbolic  conservation  laws  is  the  inherent  presence  of  discontinuities  (shocks)  in  the  solution 
(Lebacque,  1996).  Moreover,  discontinuous  weak  solutions  are  not  unique  for  hyperbolic  con¬ 
servation  laws  and  entropy  conditions  must  be  satisfied  to  obtain  physically  valid  solution 
that  is  consistent  with  human  behavior  (such  as  the  driver’s  ride  impulse)  (Ansorge,  1990; 
Velan  and  Florian,  2002).  Recently,  the  analytical  solution  for  specific  classes  of  LWR  model 
was  derived,  which  assumed  that  the  flow-density  relationship  is  governed  by  a  quadratic 
function  throughout  the  density  regime  (Wong  and  Wong,  2002b),  and  then  extended  to  the 
case  of  a  piecewise  quadratic  function  (Lu  et  al.,  2006).  Their  constructed  entropy  solutions 
are  exact  if  the  initial  condition  is  piecewise  linear  and  the  boundary  condition  is  piecewise 
constant.  The  fundamental  diagrams  in  their  works  are  continuous. 

However,  when  traffic  flow  data  are  plotted  on  the  fundamental  diagram,  the  uncongested 
and  congested  regimes  may  be  separated  by  gaps  or  discontinuities  as  shown  in  Figure  1. 
Eclie  (1961)  was  among  the  first  to  point  out  that  traffic  behaved  differently  at  different 
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density  regimes,  and  introduced  the  idea  of  a  two-regime  model  leading  to  a  discontinuous 
fundamental  diagram.  The  discontinuous  fundamental  diagrams  have  also  been  observed 
from  empirical  works  (e.g.,  Ceder,  1976;  Ceder  and  May,  1976;  Drake  et  ah,  1967;  May  and 
Keller,  1967).  In  particular,  Koshi  et  al.  (1983)  suggested  a  ’’reverse  lambda”  shape  to 
describe  the  characteristics  of  the  data  plotted  on  the  discontinuous  fundamental  diagrams. 
Further  evidence  of  the  discontinuous  fundamental  diagram  was  revealed  by  a  series  of  papers 
by  Hall  (1987),  Hall  and  Gunter  (1986)  and  Hall  et  ah  (1986).  In  addition,  Bank  (1991a, b) 
described  this  discontinuous  fundamental  diagram  as  a  two-capacity  phenomenon,  with  one 
capacity  corresponding  to  the  tip  of  the  left  leg  of  the  reverse  lambda,  and  the  other  capacity 
belonging  to  the  tip  of  the  right  leg  of  the  reverse  lambda  with  a  capacity  drop  from  the 
former  tip.  More  recently,  Cassidy  (1998)  and  Cassidy  and  Bertini  (1999)  also  confirmed 
the  capacity  drop  on  highways.  When  such  discontinuous  fundamental  diagram  is  embedded 
into  the  LWR  model,  it  is  to  our  best  knowledge  that  there  is  still  no  mathematical  theory 
on  the  existence  and  uniqueness  of  the  entropy  solutions  for  the  resultant  traffic  model. 

In  this  paper  we  assume  that  the  flow-density  relationship  q(p)  is  concave  and  is  repre¬ 
sented  by  a  piecewise  quadratic  function,  with  any  two  adjacent  pieces  joining  discontinu- 
ously  at  a  critical  density  p0.  Such  discontinuous  fundamental  diagrams  were  developed  in 
Drake  et  al.  (1967)  by  fitting  with  observed  data.  Our  procedure  to  construct  the  physically 
relevant  solutions  for  such  conservation  laws  with  discontinuous  fluxes  is  as  follows.  We  first 
explicitly  construct  the  entropy  solutions  to  a  sequence  of  approximate  problems  in  which 
the  flow-density  relationship  q(p)  is  continuous  but  tends  to  the  discontinuous  flux  when  a 
small  parameter  in  this  sequence  tends  to  zero.  We  then  explicitly  construct  the  limit  of 
the  entropy  solutions  for  this  sequence  and  consider  this  limit  solution  as  the  the  entropy 
solution  associated  with  the  discontinuous  flux.  In  order  to  verify  the  physical  relevancy 
of  such  entropy  solutions,  we  apply  our  results  to  a  few  typical  traffic  flow  examples  and 
comment  on  the  implication  of  the  solutions. 

The  organization  of  the  paper  is  as  follows.  In  Section  2  we  obtain  the  explicit  formulas 
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Figure  1:  A  typical  flow  with  two  different  concave  quadratic  functions  joining  discontinu- 
ously  at  a  critical  density  p0  with  decreasing  derivative  at  p0. 

for  the  entropy  solutions  to  Riemann  problems  with  discontinuous  flow-density  relationship, 
in  a  limit  process  involving  a  sequence  of  approximate  problems  in  which  the  flow-density 
relationship  q(p)  is  continuous  but  tends  to  the  discontinuous  flux  when  a  small  parameter  in 
this  sequence  tends  to  zero.  In  Section  3  we  obtain  the  explicit  formulas  for  the  entropy  solu¬ 
tions  with  discontinuous  flow-density  relationship  and  with  piecewise  linear  initial  condition 
and  piecewise  constant  boundary  conditions.  In  Section  4  we  provide  numerical  examples 
in  traffic  flows  to  demonstrate  the  explicit  solutions  obtained  in  Sections  2  and  3.  We  also 
compare  these  explicit  solutions  with  numerical  solutions  obtained  by  using  the  high  order 
weighted  essentially  non-oscillatory  (WENO)  schemes  (Jiang  and  Shu,  1996;  Zhang  et  ah, 
2003).  Concluding  remarks  are  given  in  Section  5. 
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2  A  sequence  of  approximate  problems  with  continu¬ 
ous  fluxes:  Riemann  problems 

The  governing  equation  for  the  LWR  model  is  the  following  scalar  hyperbolic  conservation 
law 

Pt  +  q(p)x  =  0  (1) 

with  suitable  initial  and  boundary  conditions.  Here  p  G  (0,pmax)  is  the  density,  pmax  is  the 
maximum  (jam)  density,  and  q(p )  is  the  traffic  flow  on  a  homogeneous  highway,  which  is 
assumed  to  be  a  function  of  the  density  p  only  in  the  LWR  model.  More  specifically,  the 
flow  q ,  the  density  p  and  the  equilibrium  speed  u  are  related  by 


q(p)  =  u(p)p. 


(2) 


In  this  paper,  the  flow  q(p)  is  considered  to  be  piecewise  quadratic  and  locally  concave  in 
each  piece.  Without  loss  of  generality,  we  will  concentrate  our  discussion  on  the  situation 
where  the  flow  q  is  defined  by  two  different  quadratic  functions  in  different  regimes 


where 


q(p)  = 


9i(p),  0  <  p  <  p0 

92  (p),  Po<P<  Pmax 


Flux  I:  gi  (p)  =  d0  +  dip  +  d2  p2; 


Flux  II:  g2 (p)  =  e0  +  ei  p  +  e2  p2 


(3) 

(4) 


are  two  different  quadratic  functions,  which  are  discontinuous  at  the  junction  qi(po)  >  <72(po), 
concave  in  each  piece  q"(p)  <  0  and  <72(p)  <  0,  and  concave  also  at  the  junction  q[(po)  > 
92  (Po)- 

A  typical  flow  in  this  setup  is  given  in  Figure  1.  The  general  situation  of  the  flow  q  with 
more  than  two  pieces  of  quadratic  functions  can  be  considered  with  the  same  recipe  to  each 
neighboring  pairs  of  quadratic  flow  functions. 

Since  it  is  not  known  mathematically  how  to  study  the  existence  and  uniqueness  of 
entropy  solutions  for  the  scalar  conservation  law  (1)  with  a  discontinuous  flux  function  q(p), 
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we  first  consider  the  following  sequence  of  conservation  laws  with  continuous  fluxes 


Pt  +  <f{p)x  =  0 


(5) 


where 


{qi(p)  =  d2p2  +  dip  +  d0,  0<p<po 

QmM  =  -  qi(po)~f{po+£)  ( p  -  po)  +  gi(po),  Po  <  P  <  Po  +  £ 
qUp)  =  e2 p2  +  (ei  -  2 ee2)p  +  e0  +  e2£2  -  e^,  p0  +  £  <  p  <  p„ 

with  a  simple  Riemann  type  initial  condition 


p(x,  0)  = 


Pi,  x  <  0 

Pr,  X  >  0 


Notice  that  the  narrow  £  region  connecting  the  two  pieces  of  the  discontinuous  fluxes  is 
located  to  the  right  of  the  discontinuity  point  p0.  We  could  of  course  also  construct  the 
sequence  with  the  narrow  £  region  located  to  the  left  of  the  discontinuity  point  po,  or  centrally 
around  po- 

We  remark  here  the  technical  difficulty  in  dealing  with  the  conservation  law  (1)  with 
the  flux  (6):  the  flux  is  continuous  for  £  >  0,  but  it  is  not  globally  concave.  Therefore,  we 
cannot  use  our  results  in  (Lu  et  ah,  2006)  directly.  It  is  not  possible  to  connect  the  two 
discontinuous  pieces  of  the  flux  in  a  continuous  and  globally  concave  fashion. 

We  are  interested  only  in  the  situations  that  pi  <  po  <  pr  or  pi  >  p0  >  pr,  with  the 
equality  holding  at  most  once  in  each  situation.  Otherwise,  the  Riemann  problem  would 
involve  only  one  of  the  two  quadratic  fluxes  in  (4)  and  its  solution  would  be  routine. 

Even  though  the  flux  q£(p)  is  not  globally  concave,  it  is  continuous  for  £  >  0.  Therefore, 
we  can  construct  the  solution  to  the  Riemann  problem  using  standard  techniques,  see,  e.g. 
(LeVeque,  1992),  via  convex  or  concave  hulls  from  the  graph  of  the  flux  qe(p). 


2.1  Shock:  pi  <  pr 


If  pi  <  pr,  and  without  loss  of  generality  assume  pi  <  p0  <  p0  +  £  <  pr,  then  we  construct 
the  convex  hull  of  the  set  {(p,y)  :  pi  <  p  <  pr  a nd,  y  >  q£(p)}.  Recall  that  the  convex  hull 
is  the  smallest  convex  set  containing  the  original  set.  An  illustration  of  the  convex  hull  is 
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shown  in  Figures  2  and  3  (the  shaded  region).  We  can  then  obtain  the  entropy  solution  of 
the  Riemann  problem  with  the  flux  qe(p)  as  follows,  where  £  =  x/t: 


Figure  2:  The  convex  hull  of  {(p,y)  :  pi  <  p  <  pr  and  y  >  q£(p)}  for  kAB  <  kBc- 


Let  kAB  and  kBc  be  the  slope  of  AB  and  the  slope  of  BC,  respectively,  where 


kAB  — 


g|(po  +  g)  -qi(pi) 
Po  +  e  —  Pi 


kBC  — 


Q£2(Pr)  0.2  (do  T  S) 

Pr  ~  (do  +  e) 


•  If 


kAB  <  ksc-, 


then  the  solution  is  (see  Figure  2) 


p(x,t) 


Ph  £  <  kAB 

Po  +  kAB  <  £  <  kBc 

Pr ,  £  >  kBc 


•  If 


(8) 

(9) 


kAB  >  kBC, 


(10) 
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p 


Figure  3:  The  convex  hull  of  {(p,  y)  :  pi  <  p  <  pr  and  y  >  q£(p)}  for  Pab  >  k bc- 


then  the  solution  is  (see  Figure  3) 

p(x,t)  = 


Ph 
Pr , 


C  <-  g|(Pr)-gl(ft) 

^  —  pr-pl 

C  >  <l2(Pr)-<n(pl) 
^  Pr-Pl 


(11) 


Therefore,  taking  the  limit  £  — >  0+,  we  have  the  following  explicit  solutions  to  the  Rie- 
mann  problem  with  the  discontinuous  flux  (3)  : 


•  If 


then 


If 


g2(po)  ^  qii.pl)  q2(pr)  -  q2(po) 

Po  —  Pi  Pr  —  Po 

£  <-  q2{po)~qi(pi) 

^  —  po-pi 

q2(po)-qi(pi)  <  £  <  q2(pr)-g2(po) 
PO~Pl  ,  pr-po 

£  >  q2  \Pr )  q2  \P0 ) 

^  —  Pr  —  pO 


{Ph 
Po, 
Pr, 


fe(Po)  -  <ll(pi)  >  l2(pr)  Q2  (Po) 


(12) 


(13) 


Po  -  Pf 


Pr  PO 


(14) 


then 


P  (x,t) 


(15) 


2.2  Rarefaction  wave:  pi  >  pr 

If  pi  >  pr ,  and  without  loss  of  generality  assume  pi  >  p0  +  e  >  po  >  Pr,  then  we  construct 
the  convex  hull  of  the  set  {(p,  y)  :  pr  <  p  <  pi  and  y  <  q£(p)}.  An  illustration  of  the  convex 
hull  is  shown  in  Figures  4  and  5  (the  shaded  region). 


R 


Q. 

CT 


Pr  Po  Po+E  P*  Pi 


p 


Figure  4:  The  convex  hull  of  {(p,y)  :  pr  <  p  <  pi  and  y  <  qe(p)}  for  {q^'ipi)  <  k^. 

If  we  look  at  the  upper  boundary  of  the  convex  hull  in  Figure  4,  we  can  observe  that  it 
consists  of  a  straight  line  segment  from  (p0,  q£(po))  to  (p*,  q£(pir))  (the  line  BC  in  the  figure), 
and  the  two  pieces  of  the  curve  y  =  q£(p)  for  pr  <  p  <  p0  (AB  in  the  figure)  and  p*  <  p  <  pi 
(CD  in  the  hgure),  where  BC  is  the  tangent  of  q£(p). 

For  the  slope  of  BC,  we  have  the  relationship: 
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p 

Figure  5:  The  convex  hull  of  {(p,y)  :  pr  <  p  <  pi  and  y  <  q£(p)}  for  (g|)'(p;)  >  k^. 

0.2  (P*)  9l(Po)  (  e\U  \  n  r  i  1 

-  =  (?2)  (P*)  =  2e2P*  +  ei  (16) 

P*  _  Po 

where  e2  =  e2,  ej  =  ei  —  2ee2  and  e'0  =  e0  +  e2e2  —  e\ e.  Then  we  can  get  a  quadratic  equation 
of  /?*.  The  discriminant  of  this  quadratic  equation 

A  =  +  eiPo  +  eo  —  (^2Po  +  ^iPo  +  ^o)]  >  0  (17) 

because  e2  <  0  and  e2pl  +  e[po  +  e'0  <  d2po  +  diPo  +  do  when  e  is  small  enough.  The  smaller 
root  of  (16)  is  just  what  we  want. 

Let  Pq  =  (q2)'(p*)  =  2 e'2p*  +  e\ .  Apparently,  <  (g|)'(p0)  =  e)  +  2e'2p0.  We  can  then 
obtain  the  entropy  solution  of  the  Riemann  problem  with  the  flux  qe(p )  as  follows,  where 
f  =  x/t: 
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(18) 


•  If 


as  shown  in  Figure  4,  then  the  solution  is 


pe(x,t) 


'  pu  t  f  <  (?i yi.pi) 

(viYipi)  <£<K 

<  Po,  fc§Cf<gi(po) 

9l(Po)  <  £  <  9i(Pr) 

k  Pr,  ^>g'l(Pr) 


(19) 


•  If 


(^)'(Pt)  >  k 


(20) 


as  shown  in  Figure  5,  then  the  solution  is 

( 


<&{pi)-qi(.po) 
Pi-fio 


pe(x,t)  =  { 


pi,  €  < 

SS'  CA(P0)  <  (  <  ‘A(Pr) 
l  Pi -t  C  —  (Pr) 


(21) 


Therefore,  taking  the  limit  £  — »  0+,  we  have  the  following  explicit  solutions  to  the  Ric- 
mann  problem  with  the  discontinuous  flux  (3),  where  k o  is  the  value  of  with  e  =  0: 


•  If 


then 


If 


then 


?2 (Pi)  < 


POM) 


'  Pi,  f  <  ?2(Pi) 

,  92 (Pi)  <£<ko 
<  p0,  ko  <  ^  <  q[(po) 
Q'l(po)  <  £  <  q[(pr) 
,  Pr,  £,>q'l(pr) 


q’iipi)  >  ko 


p(x,t )  =  < 


Pj, 

Po> 
g-rfl 
2d,2  ' 

Pr, 


g2(pQ~gl(po) 
Pl~P0 

Ql(Po)  <  {  <  Hl(Pr) 

?  >  9'l<Pr) 


(22) 


(23) 


(24) 


(25) 
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3  Explicit  construction  of  the  entropy  solutions 


We  now  start  the  construction  of  explicit  solutions  to  the  conservation  law  (1)  with  such 
flows  q(p),  when  the  initial  condition  is  piecewise  linear.  We  will  first  ignore  the  boundary 
conditions,  and  will  leave  the  discussion  of  the  treatment  of  piecewise  constant  boundary 
conditions  to  Sections  3.3  and  3.4.  We  begin  with  the  generalized  Riemann  problem 


P(x,  0) 


c*i  +  (3i  x,  x  <  0 
a2-\-  (32x,  x  >  0 


(26) 


We  also  assume  that,  for  the  x  range  we  are  considering,  the  initial  density  +  fax  is 
completely  contained  in  one  of  the  regimes  p  <  p0  or  p  >  p0  for  i  =  1  and  2.  This  does 
not  lose  generality,  as  we  can  break  a  single  linear  function  into  two  pieces  as  in  (26)  when 
it  crosses  the  critical  density  p0.  We  also  remark  that  we  do  not  need  to  consider  the  case 
when  both  linear  functions  cq  +  x,  for  i  —  1,2,  are  contained  in  a  single  regime  p  <  po  or 
P  >  Po,  because  this  is  covered  by  the  results  in  (Wong  and  Wong,  2002b). 

As  in  (Wong  and  Wong,  2002b  and  Lu  et  ah,  2006),  we  will  use  heavily  the  following 
simple  fact:  for  the  scalar  conservation  law  (1)  with  a  quadratic  flux  q(p)  =  a  +  b  p  +  c  p2 
and  a  linear  initial  condition  p(x,  0)  =  a  +  f3  x,  the  solution  stays  linear 


p(x,  t )  =  a(t)  +  /3(t)  x 


(27) 


with 


a  -  b(3t  (3 

aw  =  i  ,  Pi*)  = 


(28) 


1  +  2  cf3f  1  +  2  cpt 

This  simple  fact  is  the  main  reason  that  enables  us  to  obtain  explicit  formulas  for  the  entropy 
solution.  The  solution  for  each  linear  piece  of  the  initial  condition  is  given  by  (27)-(28)  until 
neighboring  waves  interact  with  each  other. 

We  now  assume  that  the  rr-axis  is  divided  into  a  number  of  elements,  within  each  of  which 
the  initial  density  is  given  by  a  linear  function  p(x,  0)  =  a  +  /3x  that  is  completely  contained 
in  one  of  the  regimes  p  <  po  or  p  >  po-  We  consider  the  solution  to  the  generalized  Riemann 
problem  with  the  two  piecewise  linear  initial  conditions.  The  left  and  right  elements  to  the 
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inner  boundary  point  under  consideration  are  denoted  by 


e=(xi,xr),  e  =  (xi,xr) 

respectively,  with  clearly  xr  =  Xi.  The  initial  condition  density  values  at  the  relevant  element 
boundaries  are  denoted  by 

Pi  =  p{xf,  0),  pr  =  p{x~i  0);  p^pix+,0),  pr  =  p(x;,0) 

We  also  denote  the  densities  in  the  left  and  right  elements  by  pi(x,  t)  and  p2(x,  t ),  respectively, 
where 

Pi(x,t)  =  ai(t)  +  fii(t)x,  p2 {x ,  f )  =  a2(t)  +  (32(t)x. 

We  discuss  the  situation  for  rarefaction  waves  in  Section  3.1  and  the  situation  for  shocks 
in  Section  3.2.  In  Sections  3.3  and  3.4  we  discuss  the  piecewise  constant  left  and  right 
boundary  conditions,  respectively.  At  last  we  provide  the  solution  procedure  in  Section  3.5. 

3.1  Case  I:  pr  >  po  >  Pi 

The  elements  e  and  e  belong  to  Flux  II  and  Flux  I  in  (4),  respectively.  Denote 

1  —  I  o  I  o  pteipo)  —  9i(Po)  /orA 

ko  —  ei  +  2e2Po  +  2e2*/ -  (29) 

V  e2 

We  have  the  following  two  sub-cases. 

3.1.1  Sub  -case  I  (a):  q'2(pr)  <  k0 

In  this  sub-case  three  new  elements  ei,  e2  and  e3  are  created  at  the  time  t  =  At,  as  shown 
in  Figure  6.  We  again  consider  only  the  time  At  smaller  than  the  smallest  time  when  the 
waves  (characteristic  lines  or  shocks)  from  the  initial  condition  intersect  with  one  another. 
The  coordinates  of  the  four  nodes  serving  as  the  end  points  of  the  three  elements  e1?  e2  and 
e3  at  time  At  can  be  determined  as 

Xi(At)  =  xr  +  q'2(pr)  At,  x2(A  t)  =  xr  +  k0At 
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x3(A  t)  =  xr  +  q[(p0)At,  x4(A  t)  =  xr  +  q'^p^At 


The  density  at  these  end  points  at  time  At  are  given  by 

p(xi(At),  At)  =  pr,  p(x2(At)~,  At)  =  ,  p(x2(At)+,  At)  =  p0, 

Ze2 

p(x3(At),  At)  =  p0,  p(x4(At),  At)  =  pt 

and  the  density  is  linear  within  each  of  the  new  elements  e,:,  i  =  1,  2,  3.  In  particular,  the 
density  within  e2  is  a  constant  p  =  po- 


Figure  6:  Sub-case  I  (a):  q2(pr)  <  ^o- 
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3.1.2  Sub-case  I  (b):  q'2(pr )  > 


In  this  sub-case,  the  situation  is  shown  in  Figure  7.  The  discontinuity  which  is  generated 
from  the  point  xr  =  Xi  will  move  to  the  right  or  left  along  a  curve.  If  the  discontinuity  moves 
to  the  left,  an  easy  way  to  determine  the  location  of  the  discontinuity  Xi  +  Ax  after  time  At 
is  through  conservation  in  the  rectangular  region  hi  with  (xi  +  Ax,  0)  and  (xi  +  q[(pl)At,  At) 
as  the  end  points  of  a  diagonal  (see  Figure  7). 


Figure  7:  Sub-case  I  (b):  q'2(pr )  >  k0. 

The  flux  at  the  left  boundary  dfli,  namely  the  number  of  vehicles  coming  from  the  left 
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boundary  into  the  region  hi  during  time  period  At  is 


rAt  rAt 

fi=  q-2\x=xr+Axdt  =  /  {e0  +  e1[a1(t)+/3i(t)(xr  +  Ax)]+e2[a1(t)+l31(t)  (xr  +  Ax)]2}dt 

Jo  Jo 

(30) 

Likewise,  the  flux  at  the  right  boundary  <9hlr,  namely  the  number  of  vehicles  leaving  the 
right  boundary  from  the  region  hi  is 


rAt 


fr  = 


Ql  \x=xi+q'1(pi)At  dt 


(31) 


cAt 


{do  +  di[a2(t)  +  j32{t){xi  +  q^p^At)}  +  d2[a2(t)  +  (J2(t){xi  +  q'l(pl)At)}2}  dt. 


The  initial  number  of  vehicles  within  the  region  hi  at  time  t  =  0  is 

px 

fb= 


r  Xr  rxi+q'^ppAt 

I  (a!i(0)  +  /5i(0)x)  dx  +  /  (ck2(0)  +  {32(0)x)  dx 

xr-\-Ax  Jxi 


and  the  final  number  of  vehicles  within  the  region  hi  at  time  t=A t  is 

^  r^l+Ql1(.Po)At  rxi+q'^tppAt 

ft=  p0dx+  (a2(At)  +  (32{At)x)  dx. 

J  xr + Ax  J  xi  H-q'J  (po )  At 

From  the  flow  conservation  principle,  we  deduce  that 


(32) 


(33) 


fl  ~  fr  +  fb  ~  ft  ~  0. 


(34) 


We  then  obtain  from  (34)  the  explicit  equation  determining  Ax  as 


Fi(Af)Aa;2  +  F2(At)Ax  +  F3(At)  =  0 


(35) 


where 

Fi(Af)  =  pi-pr 

F2(At)  =  —2[ei(pi  —  pr)At  +  2e2p0(pi  —  pr)At  +  (p0  —  pr){xi  —  xr)] 

F3(At)  =  At{el(pi  -  pr)At +  2e1pr(xr  -  xt) +  2{-2e0e2piAt  +  2die2p0ptAt 

+2d2e2plpiAt  +  2e0e2prAt  -  2d,1e2p0prAt  -  2d2e2p20prAt  -  e0xi  +  dip0xi 
+d2p20xi  -  e2p2xi  +  d0[2e2(pi  -  pr)At  +  xt  —  xr]  +  [e0  -  p0(di  +  d2p0)  +  e2p2r]xr}} 
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The  discontinuity  trajectory  can  therefore  be  determined  by  solving  the  quadratic  equation 


(35).  Then  when  F\  (At)  ^  0, 


Ax  = 


-F2(At)  -  VA 


(36) 


2F1(At) 

where  A  =  F2(At)2  —  AFi(At)F3(At).  If  the  root  Ax  is  positive,  then  the  assumption  of 
a  left-moving  discontinuity  is  incorrect.  We  can  then  repeat  the  procedure  above  assuming 
that  the  discontinuity  is  right-moving.  It  turns  out  that  the  formula  for  the  location  Ax  of 
the  discontinuity  is  still  given  by  (36).  Therefore,  we  do  not  need  to  distinguish  whether  the 
discontinuity  moves  to  the  left  or  right.  The  coordinates  of  the  three  nodes  serving  as  the 
end  points  of  the  three  elements  e3  and  e2  at  time  At  can  be  determined  as 


aq(At)  =  xr  +  Ax(At),  x2(At)  =  xr  +  q((p0)At,  x3(At)  =  xr  +  q'1(pl)At. 


The  density  at  these  end  points  at  time  At  are  given  by 


p(xl(At)+,  At)  =  p0 ,  p(x2(At),  At)  =  p0,  p(x3(At),  At)  =  pt 


and  the  density  is  linear  within  each  of  the  new  elements  e*,  i  —  1,  2.  In  particular,  the 
density  within  e!  is  a  constant  p  =  p0. 

3.2  Case  II:  pr  <  po  <  Pi 

In  this  case  the  elements  e  and  e  belonging  to  Flux  I  and  Flux  II  in  (4),  respectively.  We 
have  the  following  two  sub-cases. 

3.2.1  Sub-case  II  (a):  pr  <  p0 

(i)  If 

qi(pr)  -  g2(po)  >  92 (Pi)  -  ?2 (Po) 

Pr  -  PO  ~  Pi  ~  PO 

then  a  shock  satisfying  the  Lax  entropy  condition  (Lax,  1973)  is  generated  from  the  point 
xr  =  Xi  and  will  move  to  the  right  or  left  along  a  curve  determined  by  the  Rankine-Hugoniot 
jump  condition.  If  the  shock  moves  to  the  right,  as  shown  in  Figure  8,  an  easy  way  to 
determine  the  location  of  the  shock  xi  +  Ax  after  time  At  is  through  conservation  in  the 
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rectangular  region  f2  with  (xi,  0)  and  (xi  +  Ar,  At)  as  the  end  points  of  a  diagonal  (see  Figure 
8).  Notice  that,  since  we  consider  only  the  time  At  smaller  than  the  smallest  time  when  the 
waves  (characteristic  lines  or  shocks)  from  the  initial  condition  intersect  with  one  another, 
we  can  safely  assume  that  the  left  and  top  boundaries  of  this  rectangle,  dQi  and  <9h2t,  belong 
to  Flux  I,  and  the  right  and  bottom  boundaries  of  this  rectangle,  dQr  and  <9f2&,  belong  to 
Flux  II. 


Figure  8:  Sub-case  II  (a)  (i):  pr  <  p0  and  qi^prJ _qp2J'P0'>  >  q2^Jl2^p0'>. 


The  flux  at  the  left  boundary  dfli,  namely  the  number  of  vehicles  coming  from  the  left 
boundary  into  the  region  during  time  period  At  is 

rAt  rAt 

fi=  qi\x=xrdt=  /  [d0  +  di(ai(f)  +Pi(t)xr)  +  d2(a1(t)  +  p^t)  xr)2]  dt  (37) 

Jo  Jo 
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Likewise,  the  flux  at  the  right  boundary  <9Qr,  namely  the  number  of  vehicles  leaving  the 
right  boundary  from  the  region  12  is 


pA  t 


fr  I  Q2\x=xi+Ax  dt 

JO 
r>At 

—  /  [^o  +  ei(a;2(t)  +  (32{t){xi  +  Ax))  +  62(012(0  +  02(t)(xi  +  Arr))2]  dt-  (38) 


The  initial  number  of  vehicles  within  the  region  Q  at  time  t  =  0  is 

pxi+Ax 


ft  = 


(a2(0)  +  /32(0)x)  dx 


'X  l 


and  the  final  number  of  vehicles  within  the  region  Q  at  time  t  =  At  is 

pxr-\-Ax 

ft=  (on  (At)  +Pi(At)x)dx. 

J  xr 

From  the  flow  conservation  principle,  we  deduce  that 


(39) 


(40) 


ft  ~  fr  +  fb~ft  =  0. 


(41) 


We  obtain  the  explicit  equation  determining  Ax  by 

.Fi(Af)  Ax2  +  F2(At)  Ax  +  F3(At)  =  0  (42) 

where 

Fi(At)  =  2 A t(d2  -  e2)(pr  -  Pi)(pr  ~  Pi)  +  ixr  ~  xi)(pr  -  pt)  -  (xr  -  xt)(pr  -  pi) 

F2(At)  =  2 {-{xi  -  xr)[At(ei  +  2 e2pr)(pl  -  pr)  ~  C Pi  ~  Pr)(xt  -  ccr)}  +  diAt(pi  -  pr)[2e2At(pt  -  pr) 
+xi  -  xr\  -  2d2At(pi  -  pr)[Atei(p,  -  pr)  -  p^xi  -  xr)]} 

F3(At)  =  A t{(xi  -  xr){At[e\  +  4(d0  -  e0)e2](pl  -  pr)  +  2(d0  -  e0  -  eip;  -  e2p2)(xt  -  xr )} 

+2d2{At2[e2  +  4(d0  -  e0)e2](pi  -  pr)(p;  -  pr)  +  2At{e2\plp2r(xl  -  xr)  +  p2rpr{xr  -  xt) 
-pKpi  -  Pr)(xi  -xr)}  +  (do  -  e0  -ei Pi)(pi  -  pr)(xi  -xr)}  +  p2r(xt  -  xr)(xt  -xr)} 
-d\At{pi  -  pr)[2Ate2(p,  -  pr)  +xi-  xr\  +  2dlpr(xi  -  xr){2A te2(pl  -  pr)  +xt-  xr )} 
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The  shock  trajectory  can  therefore  be  determined  by  solving  the  quadratic  equation  (42).  If 
Fi(At)  ^  0,  then 


Ax  —  Zim^  (43) 

2F1(At)  y  J 

where  A  =  F2(At)2  —  AFi(At)  F3(At) .  If  the  root  Ax  is  negative,  then  the  assumption  of 


a  right-moving  shock  is  incorrect.  However,  as  before,  it  turns  out  that  the  formula  (43)  is 
still  valid  in  this  case. 

(ii)  If,  on  the  other  hand, 

gi(Pr)  ~  <?2(po)  q-2(pi)  ~  g2(po) 

Pr  ~  PO  Pi  ~  PO 

then  two  shocks  satisfying  the  Lax  entropy  condition  are  generated  from  the  point  xr  =  Xi 
and  one  new  element  ex  is  created  at  the  time  t  =  At,  as  shown  in  Figure  9.  We  again 
consider  only  the  time  At  smaller  than  the  smallest  time  when  the  waves  (characteristic 
lines  or  shocks)  from  the  initial  condition  intersect  with  one  another. 

By  the  Rankine-Hugoniot  jump  condition 


Ax'2{A  t)  = 


g2(a2(At)  +  (32(A t)(xt  +  Ax2(At)))  -  g2(po) 
a2(At)  +  (32(A t)(xi  +  Ax2(At))  -  p0 


If  Pi  ^  Pr,  then 


Ax2(At)  =  [e1(pl  -  pr)Aty/xr  -xi  +  2  e2p0(pl  -  pr)At^/xr  -  xt  (45) 

+  (Vxr  -  xi  -  y/2 e2(pr  -p^At  +  Xr-  Xi)(p0Xi  -  prTi  -  p0xr  +  pzxr)]/[(p;  -  pr)y/xr-xi] 


If  Pi  =  pr,  then 


Ax2(At)  =  [ei  +  e2(pi  +  p0)]At 


The  flux  at  the  left  boundary  dili,  namely  the  number  of  vehicles  coming  from  the  left 
boundary  into  the  region  during  time  period  At  is 


fi=  qi\x=xrdt=  [do  +  di(a>i(t)  +  /?i(t)  xr)  +  d2(a>i(t)  +  (3i(t)  xr)21\  dt  (47) 
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Figure  9:  Sub-case  II  (a)  (ii):  pr  <  p0  and  qi^ _q^p^  <  92 . 


The  flux  at  the  right  boundary  <9f7r,  namely  the  number  of  vehicles  leaving  the  right  boundary 
from  the  region  17  is 
/■At 

fr  /  0.2 \x=xi+Ax2  dt 

J  0 

At 

\e0  +  ei (a2(t)  +  j32(t)(xi  +  Ax2))  +  e2(a2 (t)  +  (32{t)(xi  +  Ax2))2]  dt.  (48) 


The  initial  number  of  vehicles  within  the  region  17  at  time  t  —  0  is 


fb 


fXi+Ax^ 


(a2(0)  +  (32(0)x)  dx 


'xi 


(49) 
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and  the  final  number  of  vehicles  within  the  region  12  at  time 


/•Ir+All  PXV+AX2 

ft  =  /  (ai(At)  +  /5i(At)a;)  dx  +  /  (p0)  dx 

J  xr  J  xr-\-  Axi 


(50) 


From  the  flow  conservation  principle,  we  deduce  that 


fi  ~  fr  +  fb  —  ft  —  0. 


(51) 


From  this,  we  obtain  the  explicit  equation  determining  Aaq  as 


Fi(At)  Aaq2  +  F2(At)  Ax\  +  F3(At)  =  0 


(52) 


where 


Fi(At)  =  (pr  -  pi)[2e2(pi  -pr)At+~xi-  xr] 

F2{At)  =  2 [di(pi  -  pr)At  +  2d2p0(pi  -  pr)At  +  (p0  -  pr)(xi  -  xr)][2e2(pz  -  pr)  At  +  xt  -  xr\ 

F3(At )  =  2diprAt(xi  -  xr)[2e2(pl  -  pr)At  +  xi~  xr]  +  d\{pi  -  pr)At2[2e2(pr  -  pz)Af  +  x^  -  xt\  + 
{xi  -  xr){At{e\fpl  -  pr)At  +  2e{pi(xr  -  xt)  +  2{d0[2e2(pl  -  pr)At  +  cci~  xr\ 

+e2p2(xr  -  xi)  +  e0[2e2(pr  -  Pi)At  +  xr-  Tz]}}  -  2[ei(pz  -  pr)At  +  2 e2p0(pi  -  pr)At 
+(po  -  Pi)(xi  -  xr)]Ax2(At)  +  (p,  -  pr)Ax2(At)2} 

+2d2Af{[e2  +  4(d0  -  e0)e2](pz  -  pr)(pz  -  pr)At2  +  2{e2[p,p2(a;z  -  xr)  +  p2p2(av  -  xt) 
-pKpi  -  Pr)(xt  -  xr)}  +  (do  -  eo)(Pi  -  pr)(xi  -  xr)}At  +  p2(xt  -  xr)(xi  -  xr ) 

+2(pz  -  pr)[2e2p0(pr  -  p^  At  -  (p0  -  pz)(Tz  -  xr)]Ax2(At) 

+{pi  -  Pr) (P;  -  Pr)Aa:2(At)2  -  2ei(pz  -  pr)Af [-prAx2(Af)  +  pt(xi  -xr  +  Ax2(At))]} 


Therefore,  when  F\  (At,)  0, 


Ax\ 


-F2(At)  +  VA 
2F1(At) 


(53) 


with  A  =  F2{At)2  —  4:Fi(At)F3(At). 

It  is  easy  to  find  that  the  density  within  ej  is  a  constant  p  =  p0. 

We  can  verify  again  that  the  formulas  for  the  location  Ax\  and  Ax2  of  the  two  shocks 
stay  the  same  regardless  of  whether  the  two  shocks  move  to  the  left  or  right. 
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3.2.2  Sub-case  II  (b):  pr  =  p0 

In  this  sub-case,  if  pi  =  pr,  then  one  shock  satisfying  the  Lax  entropy  condition  is  generated 
from  the  point  xr  =  Xi,  as  shown  in  Figure  8.  The  formula  for  the  location  Ax  of  the  shock 
is  the  same  as  that  given  in  (43). 

If  pi  7^  pr,  then  two  shocks  satisfying  the  Lax  entropy  condition  are  generated  from  the 
point  xr  =  Xi,  as  shown  in  Figure  9.  The  formulas  for  the  locations  Ax i  and  Ax2  of  the  two 
shocks  are  given  by  (53)  and  (45). 

3.3  Boundary  conditions  from  the  highway  entrance 

We  now  consider  the  boundary  condition  at  the  left  boundary  x  —  0,  which  is  the  highway 
entrance.  Looking  at  the  solution  (27)- (28),  we  can  see  that  the  general  solution  with  a 
linear  initial  condition  is  not  linear  in  t  for  fixed  x ,  unless  (3  =  0,  in  which  case  the  solution 
is  constant  in  t.  Therefore,  within  our  piecewise  linear  (in  space)  framework,  we  can  only 
consider  piecewise  constant  boundary  conditions. 

We  use  the  same  notation  as  used  previously.  The  interface  is  at  xr  =  Xi  =  0,  and  the  left 
density  value  pr  at  x  =  0  is  given  by  the  boundary  condition.  Otherwise,  this  is  identical  to 
the  situation  studied  in  the  two  previous  cases  for  an  internal  generalized  Riemann  problem. 
Again,  we  would  need  to  consider  the  situation  where  pr  and  the  linear  function  in  the  first 
element  with  end  values  ~pl  and  pr  belong  to  different  regimes  in  (4),  for  otherwise  the  solution 
is  the  one  obtained  in  (Wong  and  Wong,  2002b). 

3.3.1  Sub-case  B  (a):  pr  <  p0  <pl 

If  pr  <  po  <  Pi,  one  or  two  shocks  are  generated.  We  consider  only  the  situation  that  the 
shock  moves  to  the  right. 

(i)  pr  <  po 

If 

qi(pr)  -  g2(po)  >  g2(p;)  -  g2(po) 

Pr~  P0  ~  Pi  —  P0 
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a  shock  is  generated  and  its  location  Ax  after  time  At  is  given  by  (42),  except  that  the  left 
flux  fi  given  by  (37)  is  simplified  to 

fi  =  <h(Pr)  At 

and  the  top  flux  ft  given  by  (40)  is  simplified  to 

ft  =  Pr  Ax. 

Therefore,  the  coefficients  in  the  quadratic  equation  (43)  which  determines  the  shock  location 
Ax  are  simplified  to 

Fi(Af)  =  pr-p, 

F2(At)  =  2[e1(pl-pr)At  +  2e2pr(pl-pr)At+(pr-pl)(xi-xr)]  (54) 

F3(At)  =  -At{[ej  +4(d0  -  e0)e2](p,  -  pr)At  +  2(d0  ~>e0  -  -  e2pf)(xi  —  xr) 

+2dipr[2e2(pi  -  pr)At  +  xi—  xr\  +  2d2p2r[2e2(~pl  -  pr)At  +  x;  -  xr]}. 

On  the  other  hand,  if 

gl(Pr)  -  g2(p0)  q-2(pl)  -  g2(Po) 

Pr  ~  PO  Pi  ~  PO 

two  shocks  are  generated  and  their  locations  Axi  and  Ax2  after  time  At  are  given  by  (53) 
and  (45),  except  that  the  left  flux  fi  given  by  (47)  is  simplified  to 

fi  =  Qi(pr)  At 

and  the  top  flux  ft  given  by  (50)  is  simplified  to 

ft  =  Pr  Axi  +  Po  (Ax2  -  Axi). 

Therefore,  the  coefficients  in  the  quadratic  equation  (52)  which  determines  the  shock  location 
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Ax  |  are  simplified  to 


Fi  (At)  =  0 

F2(At)  =  2(pr  -  p0)[2e2(pi  -pr)At  +  xt  -xr] 

F3(At)  =  -At{2pr(dl  +  d2pr)[2e2(pi  -  pr)At  +  x;  -  xr]  +  -  pr)At 

+2e{pl(xr  -  xi)  +  2{d0[2e2(pi  -  pr)At  +  xi  -  xr]  +  e2p'f(xr  -  xj)  + 
e0[2e2(pr  -  ft) At  -  xt  +  xr}}}  +  2[ei(p,  -  pr)  At  +  2 e2p0(pi  -  pr)At 
+(po  ~  Pi)(xi  -  xr)\Ax2  +  (pr  -  pi)  Axl 

(ii)  pr  =  p0 

In  this  case  a  shock  is  generated  and  the  coefficients  in  the  quadratic  equation  (43)  which 
determines  the  shock  location  Ax  are  the  same  as  those  in  (54). 

3.3.2  Sub-case  B  (b):  pr  >  p0  >  pt 

If  pr  >  p0  >  Pi,  a  rarefaction  wave  is  formed.  We  still  use  the  notation  introduced  previously. 
The  interface  is  at  xr  =  xi  =  0,  and  the  left  density  value  pr  at  x  =  0  is  given  by  the  boundary 
condition.  Otherwise,  this  is  identical  to  the  situation  studied  in  the  two  previous  subsections 
for  an  internal  generalized  Riemann  problem.  Again,  we  would  need  to  consider  the  situation 
where  pr  and  the  linear  function  in  the  first  element  with  end  values  and  ~pr  belong  to 
different  regimes  in  (4),  for  otherwise  the  solution  is  the  one  obtained  in  (Wong  and  Wong, 
2002b). 

(i)  q'2(Pr)  A  k0.  The  formulas  are  the  same  as  those  given  in  Section  3.1. 

(ii)  q2(pr)  >  ho.  In  this  case  a  discontinuity  is  generated  and  its  location  Ax  after  time 
At  is  given  by  (35),  except  that  the  left  flux  fi  given  by  (30)  is  changed  to 

fl  =  q2<ypr)  At , 
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the  top  flux  ft  given  by  (33)  is  changed  to 

ft  =  PrAx  +  po(<?i(Po)At  -  Ax)  +  /  (a2(A t)  +  (32{A t)x)  dx, 

J  xr+q'1(po)At 

and  the  bottom  flux  ft,  given  by  (32)  is  changed  to 

^  rXr+q'tCpPAt 

fb=  (ct2(A t)  +  fi2(At)x)  dx. 

J  Xr 

Therefore,  the  coefficients  in  the  quadratic  equation  (35)  which  determines  the  shock 
location  Ax  are  changed  to 

Fi(At)  =  0 

F2(At)  =  2(po  —  Pr ) 

F3(At)  =  -2[d0  -  e0p0(di  +  d2p0)  -  pr(ei  +  e2pr)}At. 

3.4  Boundary  conditions  for  the  highway  exit 

In  this  subsection  we  discuss  the  boundary  condition  at  the  right  boundary  x  =  xenci,  which 
is  the  highway  exit. 

We  consider  a  typical  exit  setup  with  a  traffic  signal,  which  alternates  between  green  and 
red  lights.  This  is  similar  to  the  piecewise  constant  boundary  condition  considered  for  the 
entrance  in  the  previous  subsection.  The  constant  values  of  the  density  p  for  the  green  and 
red  lights  are  pz  =  0  and  p;  =  pmax,  respectively. 

When  ~Pi  =  Pmax,  corresponding  to  the  situation  of  a  red  light,  a  left-moving  shock  is 
formed.  We  again  only  consider  the  situation  pr  <  p0  <  ph  corresponding  to  the  situation 
where  the  density  at  the  left  of  xend  belongs  to  the  regime  of  Flux  I  in  (4). 

The  value  Ax  determining  the  shock  location  xend  +  Ax  after  time  At  (recall  that  in  this 
case  Ax  is  negative)  is  given  by  (43),  except  that  the  right  flux  fr  given  by  (38)  is  changed 
to 

fr  =  0, 
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the  top  flux  ft  given  by  (40)  is  changed  to 


ft  Pmax  Ax, 

and  the  bottom  flux  fb  given  by  (39)  is  changed  to 

pxr-\-Ax 


(«i(0)  +  /3i(0)x)  dx. 


Therefore,  the  coefficients  in  the  quadratic  equation  (42)  which  determines  the  shock 
location  xend  +  Ax  are  simplified  to 


Fi(At)  =  pi-pr 

F2(At)  =  -2[(h{pi  -  pr)At  +  2d2pl(pl  -  pr)At  +  (p,  -  pr)(x*  -  xr)} 

F3(At)  =  A t{d\{pi  -  pr)At  +  2d1pr(-xi  +  xr )  -  2{-[e0  +  p,(ei  +  e2pz)](xz  -  xr) 

+d0[2d2(pi  -  pr)At  +  Xi~  xr]  +  d2[-2e2p{p2lAt  +  2e2p2prAt  +  2e0(pr  -  Pi) At 
+2eip;(pr  -  pi) At  +  p2rxi  -  p2rxr}}} 

When  ~Pi  =  0,  corresponding  to  the  situation  of  a  green  light,  a  rarefaction  wave  is  formed. 
The  formulas  are  the  same  as  those  given  in  Section  3.1.1. 

We  can  also  easily  generalize  the  method  to  the  case  of  general  piecewise  constant  initial 
condition  at  the  exit,  just  as  for  the  entrance. 

3.5  Solution  procedure 

In  this  section  we  summarize  the  solution  procedure,  concentrating  on  the  discussion  of 
finding  the  earliest  time  when  the  waves  (characteristic  lines  or  shocks  or  discontinuities) 
from  the  previous  initial  condition  intersect  with  one  another  and  hence  the  construction  of 
the  entropy  solution  must  be  restarted  based  on  a  new  piecewise  linear  initial  data. 

Recall  our  assumption  that  the  x-axis  is  divided  into  a  number  of  elements,  within  each 
of  which  the  initial  density  is  given  by  a  linear  function  p(x,  0)  =  a  +  f3x  that  is  completely 
contained  in  one  of  the  regimes  p  <  po  or  p  >  p0.  We  consider  several  cases  of  wave 
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interactions  separately  below,  and  then  take  the  smallest  time  from  these  cases,  which  will 
serve  as  the  time  to  restart  the  solution  procedure  with  a  new  piecewise  linear  initial  data. 

Once  the  smallest  time  of  wave  interactions  is  found,  the  density  function  is  obtained  at 
this  time  using  the  previous  formulas  as  a  new  piecewise  linear  function,  and  the  procedure 
is  repeated,  until  the  final  desired  time  is  reached. 

For  an  element  e  in  which  the  initial  linear  density  profile  is  an  increasing  function,  the 
natural  break  time,  can  be  determined  by  the  following  formula  (Whitham,  1974): 

Te  =  ~^p)i 

Notice  that  under  our  assumption  (concave  quadratic  flux  q  and  increasing  linear  density 
p)  the  denomination  is  a  negative  constant,  hence  re  is  a  positive  constant.  Indeed,  if  the 
element  e  =  ( xi,xr )  and  the  initial  condition  density  values  at  the  element  boundaries  are 
pi  =  p{x+ , 0)  and  pr  =  p(x~ ,  0),  and  assuming  that  the  flux  function  in  the  element  e  is: 

q(p)  =  a0  +  ai  p  +  a2  p2 


then 

Te  =  /  A-  >  Pl  (55) 

(  oo  otherwise 

We  now  consider  the  situation  that  two  adjacent  nodes  xm  and  xm+\  correspond  to  a  pair 
of  adjacent  characteristics  and/or  discontinuities  which  will  intersect  at  time  At.  This  is  the 
most  difficult  case  since  no  closed  form  formula  exists  for  the  intersecting  time  At.  We  will 
therefore  resort  to  a  nonlinear  equation  solver  such  as  the  Newton’s  method. 

We  denote  Gm(At )  =  Axm.  Likewise,  the  location  of  the  trajectory  from  the  node  xm+i 
at  time  At  is  xm+\  +  Axm+\  and  we  denote  Gm+i(Af)  =  Arm+i.  The  displacements  Gm(At) 
and  Gm+i  (At)  are  determined  by  (43)  or  (53)  or  (46).  Therefore,  we  can  define  the  function 

S (At)  =  xm+i  -xm  +  Gm+i(At)  -  Gm(At)  (56) 


which  measures  the  distance  between  the  two  shocks.  Clearly,  S'(O)  =  xm+\  —  xm  >  0,  and 
we  would  like  to  find  the  root  of  S(At),  which  corresponds  to  the  time  that  the  two  shocks 


intersect.  As  for  each  fixed  At,  S(At)  and  S' (At)  can  both  be  readily  computed,  we  can  easily 
set  up  a  Newton  iteration  to  solve  for  the  root  of  S(At)  with  At  =  0  as  the  initial  guess.  We 
can  follow  the  discussion  in  the  appendix  of  Lu  et  al.  (2006)  to  discuss  the  uniqueness  of  the 
solution  to  equation  (56)  before  the  natural  break  time,  thereby  facilitating  the  convergence 
of  the  Newton  iteration  process. 


4  Numerical  examples 


In  this  section  we  provide  three  numerical  examples  to  illustrate  the  explicit  formulas  for  the 
entropy  solutions  obtained  in  the  previous  sections.  The  flow-density  relationship  is  given 
by: 


/  \  __  f  qi(p)  =  -0.4p2  +  lOOp,  0  <  p  <  50 

4[f  )  \  q2(p)  =  -0.02p2  -4p  +  3850,  50  <  p  <  350. 

see  Figure  1.  The  units  of  flow  and  density  are  expressed  in  veh/h  and  veh/km,  respectively. 

The  discontinuity  in  the  fundamental  diagram  occurs  at  p  =  50  veh/km,  and  the  flows  on 

the  immediate  left  and  right  of  this  continuity  are  q  =  4000  veh/h  and  q  =  3600  veh/h, 

respectively,  with  a  capacity  drop  of  400  veh/h. 

We  also  use  the  fifth  order  finite  difference  WENO  scheme  (Jiang  and  Shu,  1996;  Lu  et 
al.,  2006;  Zhang  et  ah,  2003)  to  compute  the  solution  and  make  a  comparison.  The  purpose 
of  this  comparison  is  two  fold:  first  to  validate  the  computation  of  WENO  schemes  against 
the  presumably  exact  entropy  solutions  obtained  by  the  procedure  in  this  paper;  and  second 
to  demonstrate  that  high  resolution  numerical  schemes  such  as  the  WENO  schemes  with 
adequate  numerical  viscosities  do  converge  to  the  same  entropy  solutions  that  we  obtain  an¬ 
alytically  in  this  paper  for  discontinuous  fluxes.  We  remark  that  the  choice  of  the  viscosity 
coefficient  in  the  Lax-Friedrichs  building  block  of  the  WENO  scheme  must  be  chosen  pro¬ 
portionally  to  -A  when  the  initial  condition  crosses  the  discontinuity  of  the  flux  p  =  Po  =  50, 
which  is  consistent  with  the  analytical  procedure  adopted  in  this  paper  to  derive  entropy 
solutions,  namely  to  use  a  sequence  of  approximate  solutions  with  continuous  fluxes  with  pro¬ 
gressively  sharper  gradients  and  to  take  its  limit  to  define  the  entropy  solution  corresponding 
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to  the  discontinuous  flux. 


The  first  and  second  numerical  examples  are  generalized  Riemann  problems  with  the 
following  initial  and  boundary  conditions 


pO,o) 


+  flix,  0  <  x  <  10 
a2  +  P2X,  10  <  x  <  20 


p(0,t)  =  p(0, 0). 


4.1  Example  1  (Shocks) 


Example  l(a-l):  =  20,  /?i  =  0,  a?2  =  200,  /?2  =  0 


We  have 


Example  l(a-2):  a\ 


g2(50)-gi(20)  g2(200)-g2(50) 

50  -  20  200  -  50 

15,/?!  =  0.5,  a2  =  200,  fa  =  1 


We  have 

g2(50)-gi(20)  g2(210)-g2(50) 

50  -20  210  -50 

The  solutions  are  shown  in  Figure  10. 


Figure  10:  The  exact  entropy  solution  obtained  by  the  procedure  in  Section  3  (solid  line) 
and  the  numerical  solution  obtained  by  using  WENO  scheme  with  N  =  1000  uniform  grid 
points  (circles)  at  the  time  t=30  min  for  Example  1  (a-1)  and  (a-2),  respectively. 


Example  1  (b-1):  op  =  45,/?!  =  0,  a2  =  80,  fJ2  =  0 
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We  have 


g2(50)  —  <?i  (45)  g2(80)-g2(50) 

50  -  45  80  -  50 

Example  1  (b-2):  ay  =  50,  /?i  =  — 0.5,a2  =  80,  /32  =  —1 
We  have 

g2(50)-g1(45)  g2(70)-g2(50) 

50  -45  <  70  -  50 

The  solutions  are  shown  in  Figure  11. 


Figure  11:  The  exact  entropy  solution  obtained  by  the  procedure  in  Section  3  (solid  line) 
and  the  numerical  solution  obtained  by  using  WENO  scheme  with  N  =  1000  uniform  grid 
points  (circles)  at  the  time  t=12  min  for  Example  1  (b-1)  and  (b-2),  respectively. 


The  exact  solution  of  these  problems  can  be  worked  out  using  the  procedure  in  this 
paper,  shown  as  solid  lines  in  Figures  10  and  11,  in  comparison  with  the  numerical  solution 
obtained  by  the  WENO  scheme  using  N  =  1000  uniform  grid  points,  shown  as  circles.  In 
both  cases,  we  can  see  that  the  two  results  agree  very  well. 

4.2  Example  2  (Rarefaction  waves) 

In  this  case  ko  =  —11.6569  in  (22)-(25). 

Example  2  (a-1):  =  300,  =  0,  a2  =  20,  /32  =  0 
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<72  (300)  <  k0  and  the  solution  is  shown  in  Figure  12. 
Example  2  (a-2):  cti  =  300,  f3\  —  1,  a.2  —  30,  fa  =  —1 
92(310)  <  k0  and  the  solution  is  also  shown  in  Figure  12. 


Figure  12:  The  exact  entropy  solution  obtained  by  the  procedure  in  Section  3  (solid  line) 
and  the  numerical  solution  obtained  by  using  WENO  scheme  with  N  =  1000  uniform  grid 
points  (circles)  at  the  time  t=6  min  for  Example  2  (a-1)  and  (a-2),  respectively. 


Example  2  (b-1):  ot\  =  80,  (3\  =  0,  «2  =  20,  fh  —  0 
<72(80)  >  k0  and  the  solution  is  shown  in  Figure  13. 

Example  2  (b-2):  a±  =  80,  (3\  =  — 1,  «2  =  30,  /?2  =  —1 
<72(70)  >  k0  and  the  solution  is  also  shown  in  Figure  13. 

The  exact  solution  of  these  problems  can  be  worked  out  using  the  procedure  in  this  paper, 
shown  as  solid  lines  in  Figures  12  and  13,  in  comparison  with  the  numerical  solution  obtained 
by  the  WENO  scheme  using  N  =  1000  uniform  grid  points,  shown  as  circles.  Again,  in  both 
cases,  we  can  see  that  the  two  results  agree  very  well. 
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Figure  13:  The  exact  entropy  solution  obtained  by  the  procedure  in  Section  3  (solid  line) 
and  the  numerical  solution  obtained  by  using  WENO  scheme  with  N  =  1000  uniform  grid 
points  (circles)  at  the  time  t=6  min  for  Example  2  (b-1)  and  (b-2),  respectively. 

4.3  Example  3  (Wide  moving  jam) 

Consider  a  long  homogeneous  freeway  of  length  30  km.  We  now  assume  the  following  initial 
density 

,  ,  _  f  250,  15  <  x  <  17 

^  '  (  pd,  otherwise 

which  represents  the  traffic  condition  after  an  incident  (recurrent  or  non-recurrent)  at  x  —  17. 
The  left-hand  entrance  of  the  highway  is  always  kept  at  a  density  of  pd  all  time.  We  consider 
three  cases:  (a)  pd  =  50,  (b)  pd  =  45,  and  (c)  pd  =  40  veh/km,  with  entrance  (demand)  flows 
of  4000,  3690,  and  3360  veh/h,  respectively.  The  demand  of  the  first  two  cases  exceed  the 
capacity  of  Flux  11  (3600  veh/h),  whereas  the  demand  of  the  third  case  is  below  this  capacity. 
Figures  14-16  plot  the  exact  entropy  solution  (solid  line)  and  the  numerical  solution  obtained 
by  the  fifth  order  WENO  finite  difference  scheme  with  N  =  1000  uniform  grid  points  (circles) 
for  all  three  cases.  We  can  see  that  the  two  results  agree  very  well. 

From  the  figures,  we  can  clearly  see  the  formation  of  a  wide  moving  jam  with  two  very 
sharp  shock  fronts  on  both  ends  of  the  moving  jam,  which  is  a  traffic  phenomenon  that  is 
commonly  observed  on  highways  (Kerner  and  Rehborn,  1996),  and  was  well  studied  using 
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the  three-phase  traffic  theory  (Kerner,  2004),  cell  transmission  model  (Lin  and  Lo,  2003), 
and  traveling  wave  analysis  (Jin  and  Zhang,  2003;  Kerner  and  Konhauser,  1994;  Zhang  and 
Wong,  2006;  Zhang  et  al.,  2006).  In  Examples  3  (a)  and  3  (b),  when  the  demand  from  the 
left-hand  highway  entrance  is  higher  than  the  capacity  of  the  congested  regime  (Flnx  II  in 
this  example),  the  wide  moving  jam  lasts  for  more  than  30  minutes,  whereas  in  Example  3 
(c)  when  the  demand  is  lower  than  this  reduced  capacity,  the  jam  dissolves  very  quickly.  It 
is  interesting  that  the  wide  moving  jam  can  also  be  formed  using  the  first  order  LWR  model, 
and  does  not  require  a  linear  flow-density  relationship  in  the  congested  regime  as  discussed 
in  Lin  and  Lo  (2003). 

We  also  show  in  Figure  17  the  time-space  diagrams  and  three-dimensional  speed  plots  of 
these  three  cases  to  illustrate  the  results. 

5  Conclusions 

In  this  paper  we  consider  the  explicit  construction  of  physically  relevant  entropy  solutions  of 
a  class  of  conservation  laws  with  discontinuous  flux  functions,  which  are  piecewise  quadratic 
and  locally  concave  in  each  piece.  We  treat  this  problem  as  the  limit  of  a  sequence  of  approx¬ 
imate  problems  in  which  the  fluxes  are  continuous  functions  but  with  progressively  sharper 
gradients.  We  have  presented  explicit  formulas  for  such  entropy  solutions  for  both  the  simple 
Riemann  initial  conditions  and  for  piecewise  linear  initial  conditions  and  piecewise  constant 
boundary  conditions.  We  demonstrate  these  explicitly  constructed  entropy  solutions  to  rep¬ 
resentative  traffic  flow  problems  and  compare  them  with  numerical  solutions  obtained  with 
high  order  WENO  schemes. 
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Figure  14:  The  exact  entropy  solution  (solid  line)  and  the  numerical  solution  obtained  by 
WENO  scheme  (circles)  for  Example  3  (a)  pd  =  50.  Top  left:  t  —  5  min;  top  right:  t  —  15 
min;  bottom  left:  t  —  30  min;  bottom  right:  t  =  60  min. 
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Figure  15:  The  exact  entropy  solution  (solid  line)  and  the  numerical  solution  obtained  by 
WENO  scheme  (circles)  for  Example  3  (b)  p ^  =  45.  Top  left:  t  —  5  min;  top  right:  t  —  15 
min;  bottom  left:  t  —  30  min;  bottom  right:  t  =  60  min. 
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Figure  16:  The  exact  entropy  solution  (solid  line)  and  the  numerical  solution  obtained  by 
WENO  scheme  (circles)  for  Example  3  (c)  pd  =  40.  Top  left:  t  —  5  min;  top  right:  t  —  15 
min;  bottom  left:  t  —  30  min;  bottom  right:  t  —  60  min. 
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Figure  17:  The  exact  entropy  solutions  for  Example  3.  Top:  case  (a);  middle:  case  (b); 
bottom:  case  (c).  Left:  time-space  diagram;  right:  3D  speed  plot. 
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